NUMERICAL  ANALYSIS  OF  FILTER  BASED  STABILIZATION  FOR  EVOLUTION 

EQUATIONS 

VINCENT  J.  ERVIN*,  WILLIAM  J.  LAYTONt,  AND  MONIKA  NEDA* 

Abstract.  We  consider  filter  based  stabilization  for  evolution  equations  (in  general)  and  for  the  Navier-Stokes  equations 
(in  particular).  The  first  method  we  consider  is  to  advance  in  time  one  time  step  by  a  given  method  and  then  to  apply 
an  (uncoupled  and  modular)  filter  to  get  the  approximation  at  the  new  time  level.  This  filter  based  stabilization,  although 
algorithmically  appealing,  is  viewed  in  the  literature  as  introducing  far  too  much  numerical  dissipation  to  achieve  a  quality 
approximate  solution.  We  show  that  this  is  indeed  the  case.  We  then  consider  a  modification:  Evolve  one  time  step,  Filter, 
Deconvolve  then  Relax  to  get  the  approximation  at  the  new  time  step.  We  give  a  precise  analysis  of  the  numerical  diffusion 
and  error  in  this  process  and  show  it  has  great  promise,  confirmed  in  several  numerical  experiments. 


1.  Introduction.  Simulations  in  critical  settings  often  struggle  with  numerical  artifacts  created  by 
under  resolution  of  the  spacial  scales  and  temporal  dynamics  in  the  model,  e.g.,  Brown  and  Minon  [BM95]. 
Often,  as  soon  as  there  are  sufficient  computational  resources  for  full  resolution,  the  demands  of  the  appli¬ 
cation  force  coupling  to  other  processes,  making  the  target  simulation  again  under  resolved.  The  needs  of 
under  resolved  simulations  have  led  to  a  number  of  stabilizations  in  Computational  Fluid  Dynamics.  Herein 
we  consider  a  filter  based  stabilization  extending  the  idea  of  “  evolve  one  time  step  then  filter” .  This  idea 
was  developed  by  Boyd  [Boyd98],  Fischer  and  Mullen  [FM01],  [MF98],  and  used  by  Dunca  [D02],  Mathew  et 
al.  [MLF03]  made  the  important  connection  that  “evolve  then  filter”  induces  a  new  implicit  time  relaxation 
term  into  the  discretization  that  acts  to  damp  oscillations  in  marginally  resolved  scales.  The  connection  to 
time  relaxation  links  it  to  work  of  Schochet  and  Tadmor  [ST92],  Roseneau  [R89],  Adams,  Kleiser,  Leonard 
and  Stolz  [AL99],  [AS01],  [AS02],  [SA99],  [SAKOla],  [SAKOlb],  [SAK02],  Dunca  [D02],  [D04],  [DE04]  and 
[LN07],  [LMNR06] ,  [ELN07],  [L07b],  [CL08]. 

To  present  this  connection,  consider  the  explicit  method  for  a  nonlinear  evolution  equation 


du 

dt 


+  F(u )  =  0  . 


Let  overbar  denote  a  local  spacial  filter  with  a  filter  radius  (ultimately  related  to  the  spacial  mesh  width) 
S.  Thus,  u  denotes  the  “well  resolved”  part  of  u  and  u!  =  u  —  u  denotes  the  marginally  resolved  part  of  u. 
Add  to  the  explicit  method  one  uncoupled  filter  step.  The  method  we  extend,  and  then  analyze,  adds  one 
uncoupled  filter  step  to  stabilize  an  explicit  method:  given  un  compute  un+l  by 

wn+l  —  nn 

Stepl:  Compute  wn+1  via:  - — - b  F(un)  =  0,  (1.1) 

Step  2:  Filter  wn+1  to  obtain  un+1 :  un+1  =  (1.2) 


Both  steps  can  be  done  by  black  box  modules.  The  consistency  error  of  (1.1),  (1-2)  (the  error  made 
after  one  time  step,  starting  with  exact  values)  is  0(At2  +  S 2)  provided  the  filter  is  second  order  (i.e., 
u  =  u+0(52)).  This  suggests  a  global  error  of  0(At  +  Jy ).  Following  Mathew  et  al.  [MLF03],  eliminating 
Step  2  gives 


,n+ 1  _  yU 

At 


+  F{un)  +  —  (w;n+1  -  wn+l)  =  0  , 


(1.3) 


which  is  a  time  relaxation  discretization  of  the  original  problem  with  time  relaxation  coefficient  1/A t. 

Suppose  in  addition  that  the  filter  is  a  differential  filter,  Germano  [Ger86],  i.e.,  that  w(x,  t)  is  the  solution 
of 


—S2  Aw  +  w  =  w. 


^Department  of  Mathematical  Sciences,  Clemson  University,  Clemson,  SC,  29634-0975,  USA.  email:  vjervin@clemson.edu. 
Partially  supported  by  the  US  Army  Research  Office  under  grant  W911NF-05- 1-0380. 

^Department  of  Mathematics,  University  of  Pittsburgh,  Pittsburgh,  PA,  15260,  USA.  email:  wjl@pitt.edu.  Partially 
supported  by  the  NSF  under  grant  DMS-0810385. 

•^Department  of  Mathematical  Sciences,  University  of  Nevada  Las  Vegas,  Las  Vegas,  NV,  89154-4020,  USA.  email: 
monika.neda@unlv.edu  . 


1 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

JAN  2010  2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2010  to  00-00-2010 

4.  TITLE  AND  SUBTITLE 

Numerical  Anlaysis  of  Filter  Based  Stabilization  for  Evolution  Equations 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Pittsburgh, Learning  Research  and  Development 

Center ,3939  O’Hara  Street, Pittsburgh, PA, 15260 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

see  report 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF:  17.  LIMITATION  OF 

_ _ _  ABSTRACT 

18.  NUMBER  19a.  NAME  OF 

OF  PAGES  RESPONSIBLE  PERSON 

a.  REPORT  b.  ABSTRACT  c.  THIS  PAGE  Same  US 

unclassified  unclassified  unclassified  Report  (SAR) 

24 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


It  then  follows  from  w  —  w  =  —52Aw,  and  un+1  =  wn+1,  that  (1.1), (1.2)  becomes 


un+l  _  un 

At 


F(un) 


s2 

At 


A un+l  =  0 , 


which  is  an  Implicit-Explicit  time  discretization  of  the  artificial  viscosity  method  [ARW95].  The  artificial 
viscosity  coefficient  is  S2 /At  which  can  result  in  low  accuracy  and  large  amounts  of  numerical  diffusion 
depending  on  the  relative  scalings  of  At  and  6. 

When  used  with  implicit  methods,  the  effects  of  filter  based  stabilization  are  less  clear.  The  goals  of  this 
report  are  (i)  to  understand  their  effect  when  used  with  implicit  time  stepping  methods,  (ii)  to  increase  the 
accuracy  of  filter  based  stabilization,  and  (iii)  to  decrease  the  large  amounts  of  numerical  diffusion  implicitly 
induced  by  the  filter  step.  We  characterize  the  numerical  diffusion  induced  by  the  various  combinations  of 
filter,  deconvolution,  relaxation  stabilization  beginning  in  Theorem  1.2.  In  Sections  3  through  5  we  apply 
the  stabilization  to  the  Navier-Stokes  equations  and  show  how  the  induced  numerical  diffusion  is  reduced 
by  relaxation.  We  give  a  complete  stability  and  error  analysis  for  these  stabilizations  for  the  Navier-Stokes 
equations. 

1.1.  Filter  based  stabilization  of  implicit  methods.  The  fundamental  problem  of  filter  based 
stabilization  is  that,  if  not  carefully  done,  it  reduces  accuracy  by  introducing  large  amounts  of  numerical 
diffusion.  We  consider  next  the  more  complex  case  of  implicit  methods  with  filtering  (to  stabilize),  decon¬ 
volution  (to  increase  accuracy)  and  relaxation  (to  reduce  numerical  diffusion).  Let  V  <— >  L  V'  be  Hilbert 
spaces  with  duality  pairing  an  extension  of  the  L-inner  product,  denoted  by  <  ■,  •  >l,  and  j|  •  ||l  its  associated 
norm.  Let  F  :  V  — >  V'  satisfy 


<  F{v),v  >l=  0  ,  for  all  v  €  V. 

Consider  the  exactly  conservative  evolution  equation:  Find  u  :  [0,T]  — +  V  satisfying  u(0)  =  u°  £  V  and 

du 

—  +  F(u)  =  0  for  t  >  0. 

Since  this  equation  is  exactly  conservative  under  the  above  condition  on  F,  any  dissipation  of  energy  is  a 
numerical  artifact.  Let  the  filter  G  :  V'  — >  V,  and  deconvolution  operator  D  :  V  —>  V,  be  bounded,  self 
adjoint  and  positive  linear  operators,  i.e.  SPD,  which  satisfy  the  minimal  compatibility  conditions  that  as 
operators  :  L  — >  L 


DG  and  (/  —  DG)  are  L  self  adjoint  and  positive. 


The  van  Cittert  deconvolution  operators  defined  in  Section  2.3  and  many  others  satisfy  (1.4). 

Algorithm  1.1.  Pick  %  €  [0,1]  and  At  >  0. 

Step  1:  Given  un  find  wn+1  satisfying 


At 


F(- 


n+l 


-)  =  o. 


Step  2:  (a)  Filter:  Compute  wn+1  =  G(wn+1) 

(b)  Deconvolve:  Compute  D(wn+1) 

(c)  Relax: 

un+ 1  :=  (1  -  x)wn+l  +  xD(/uA+i) 


(1.4) 


(1.5) 


Relaxation  was  introduced  into  filter  based  stabilization  by  Fischer  and  Mullen  [FM01],  [MF98]  to  reduce 
the  induced  numerical  diffusion.  We  introduce  deconvolution  in  Step  2  to  increase  accuracy,  see  Theorem 
5.4  in  Section  5,  and  compare  Table  6.1  (van  Cittert  deconvolution,  see  Definition  2.2)  with  Table  6.2  (no 
deconvolution);  see  also  Table  6.4.  We  quantify  the  numerical  diffusion  introduced  by  Step  2  next. 


Theorem  1.2.  Let  u  =  (1  —  x)w  +  xD{w).  Then 

I H  li  -  I M  ll  =  x(2  -  X)  {{I  -  DG)w,  w)L  +  x2  (( I  ~  DG)w,  DGw ) L 


2 


(1.6) 


An  approximate  solution  given  by  Algorithm  1.1  satisfies,  for  any  l  >  0, 

I  Ml  ll  +  At  -  DG>"’  w")l  +  £  ((!  ~  DG)wn,  DGwn)L 

n— 0  ^ 

The  numerical  dissipation  introduced  by  Step  2  at  each  time  step  is 

X(2a~X)  (0  -  DG)wn,  wn)L  +  ^  ((/  -  DG)wn ,  DGwn) L  >  0.  (1.7) 

Proof.  Take  the  L  inner  product  of  Step  1  with  (ierl+1  +  un ),  use  <  F(v),v  >l  =  0,  and  rearrange  the 
result.  This  gives 

^  (IK+1ll!  -  \\un\\l)  +  T  [lkn+1|||  -  IK+1|||]  =  o. 

Step  2c  can  be  rearranged  to  read 

u  +  x(I  ~  DG)w  =  w  . 

Take  the  L  inner  product  with  w.  This  gives 

(u,w)L  +  x((I  ~  DG)w,w)l  =  \\w\\2l. 

The  second  term  is  nonnegative  x  {(I  ~  DG)w,  w) L  >  0.  Apply  the  polarization  identity  to  the  first  term. 
This  gives,  after  simplification 

I  Mil  -  \  W-  w\\l  +  2y((J-  DG)w,w)l  =  IMli- 

Step  2c  can  also  be  rearranged  to  read 

u  —  w  =  — x(7  —  DG)w ,  so  ||u  —  w\\\  =  y2  ((/  —  DG)w,  (I  —  DG)w)L  . 

Thus, 

IMI  ll  =  I  lMl  ll  +  x(2  -  x)  {(!  -  DG)w,  w)l  +  x2  {(I  -  DG)w,  DGw) L  , 

which  is  the  first  claim.  Inserting  this  with  u  =  un+1,w  =  wn+1  gives 

^  (IK+1||!  -  ||«n|ll)  +  ^  [x(2  -  x)  ((I  -  DG)wn+1 ,  wn+1) L  +  x2  ((I  DG)wn+1 ,  DGwn+1) L]  =  0. 

Summing  gives  the  energy  estimate.  That  the  term  in  brackets  is  nonnegative  and  thus  dissipated  energy 
follows  since  (/  —  DG)  is  SPD  for  the  first  term.  The  second  term  is  nonnegative  since  DG  and  /  —  DG 
commute  and  are  both  SPD. 

□ 

The  form  of  the  induced  numerical  diffusion  suggests  the  natural  scaling  (used  in  our  tests  in  Section  5) 
of  the  relaxation  parameter 


X  ~  O(At). 

With  this  scaling  the  second  term  in  the  expression  for  the  numerical  dissipation  is  a  higher  order  term;  the 
numerical  dissipation  is  dominated  by  the  first  term.  On  the  other  hand,  fluid  flow  problems  are  complex 
and  the  above  suggestive  scaling  is  not  likely  to  be  optimal  in  general.  The  determination  of  parameter 
choice  better  than  the  above  (and  which  can  possibly  be  applied  locally)  is  an  open  problem.  For  a  properly 
chosen  deconvolution  operator,  (/  —  DG)w  can  be  very  small  for  smooth  w,  so  deconvolution  reduces  the 
numerical  diffusion  added  by  filtering  to  the  large  scales. 
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Fig.  1.1.  ( I  —  DG)  for  N  =  0:  thin  line,  N  =  1:  thick  line,  and  DG(I-DG):  broken  line. 


Remark  1.3  (Interpreting  the  numerical  dissipation  term  in  Theorem  1.2).  For  the  periodic  problem 
L2(0,2tt),  G  a  differential  filter  (  G  =  (— <52 A  +  l)^1  under  periodic  boundary  conditions),  D  a  van  Cit- 
tert  deconvolution  operator,  Definition  2.2,  the  action  of  the  numerical  dissipation  in  Theorem  1.2  can  be 
calculated  wavenumber  by  wavenumber.  Let  w(x)  =  w(k)e~lkx .  The  first  two  van  Cittert  deconvolution 
operators  are  Do  =  I,  and  D\  =  21  —  G.  For  the  first  term  in  the  numerical  dissipation  (1.7)  we  then  have 

(( I -DG)w,w )  =  Y, 

k 


We  calculate 


(I-D0G)(k) 

{I^DlGffk) 


(6k)2 
(5k)2  +  1’ 

(  m2  y 

U6fc)2  +  iy 


We  plot  these  in  Figure  1.1  (along  with  the  analogous  graph  for  the  second  term  in  the  numerical  diffusion 
expression  (1.7)  for  N  =  0,  DqG  (I  —  DoG)(k)  =  (6k)2 /  ((6k)2  +  1 )  ).  Note  that  the  numerical  diffusion 
increases  as  the  wavenumber  (length  scale)  increases  (deceases).  Increasing  the  order  of  deconvolution  (here 
only  from  N  =  0  (thin  line)  to  N  =  1  (thick  line))  enhances  this  effect.  Note  also  that  the  dotted  line 
(representing  the  second  term)  is  considerably  smaller  than  the  other  terms,  before  even  considering  its 
smaller  coefficient. 

2.  Preliminaries.  We  consider  the  above  algorithms  for  the  Navier-Stokes  equations  (NSE)  in  a  poly¬ 
hedral  domain  £2  : 


(2.1) 
(2.2) 

Two  boundary  conditions  are  considered:  either  no  slip 

u  =  0  on  dtt  for  t  >  0, 


ut  +  u  ■  Vu  —  uAu  +  Vp  =  /  ,  in  fl  C  (d  =  2, 3),  t  >  0, 
V  •  u  =  0  ,  in  fl,  t  >  0, 
u(a;,0)  =  u°(x),  in  £2. 
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(2.3) 


or  under  periodic  with  zero-mean  boundary  conditions.  In  this  later  case  f l  =  (0 ,L)d  and 

u(x  +  Lej,t)  =  u(x,t)  j  =  l,-’*-,d  and, 

/  (f>dx  =  0  for  (f>  =  u,  u°,  /,  and  p. 

Jn 

There  are  differences  between  filtering  under  no  slip  or  periodic  boundary  conditions.  Under  periodic 
conditions,  incompressibility  is  preserved  by  the  simple  differential  filter: 


—52Aw  +  w  =  w. 


(2.4) 


Under  no  slip  boundary  conditions  incompressibility  is  not  preserved  by  the  above  simple  differential  filter. 
It  must  be  replaced  by  the  Stokes  differential  filter.  Given  w(x,t),  w(x,t)  is  the  solution  of 


— 6" Aw  +  w  +  VA  =  w ,  and  V  •  w  =  0  in  Q, , 
W  =  0  on  dtt. 


(2.5) 


2.1.  Notation.  Until  this  point  we  have  suppressed  the  (secondary)  spacial  discretization.  To  go 
further  it  must  be  specified  with  accompanying  technical  points.  The  L2(U)  norm  and  inner  product  will  be 
denoted  by  ||-||  and  (•,•).  Likewise,  the  LP(U)  norms  and  the  Sobolev  W*(Q)  norms  are  denoted  by  ||  •  \\lp 
and  ||  •  || Wk,  respectively.  For  the  semi-norm  in  Wk(Q)  we  use  |  •  \Wk.  Hk  is  used  to  represent  the  Sobolev 
space  Wj  (9),  and  ||  •  ||fc  denotes  the  norm  in  Hk .  The  space  H~k  denotes  the  dual  space  of  Hq.  For  functions 
v(x,i)  defined  on  the  entire  time  interval  (0,T),  we  define  (1  <  m  <  oo) 


IMloo ,fc  :=  EssSup[0'T]\\v(t, -)|U  ,  and  ||u| 


m,k  •  — 


1/m 


v(t,-Wdt 


We  base  our  analysis  on  the  finite  element  method  (FEM)  for  the  spacial  discretization  (and  believe  that 
the  results  extend  to  many  other  variational  methods) .  To  begin,  if  we  are  in  the  case  of  periodic  with  zero 
mean  boundary  conditions,  let  the  velocity  space  X  be  the  periodic  H 1  vector  functions  with  zero  mean 
and  the  pressure  space  Q  the  L2  functions  with  zero  mean: 

X  :=  n  L20(n))d,  Q  :=  L2(U). 

In  the  case  of  no  slip  boundary  conditions  the  analogous  choice  is 

X~  \md,  Q:=L2(n). 

We  use  as  the  norm  on  X  the  H 1  semi-norm  which,  because  of  the  boundary  condition,  is  a  norm,  i.e. 
for  v  €  X,  ||d||  x  :=  l|Vu||L2.  We  denote  the  dual  space  of  X  by  X*,  with  the  norm  ||  •  ||*.  The  space  of 
divergence  free  functions  is  given  by 

V:={vGX  :  (V-v,q)  =  0  Vg  G  Q}  . 

We  shall  denote  conforming  velocity,  pressure  finite  element  spaces  based  on  an  edge  to  edge  triangula¬ 
tions  of  fi  (with  maximum  triangle  diameter  h)  by 

Xh  CX,  QhC  Q. 

We  shall  assume  that  Xh,  Qh  satisfy  the  usual  inf-sup  condition  necessary  for  the  stability  of  the  pressure, 
e.g.  [G89] .  The  discretely  divergence  free  subspace  of  Xh  is 

Vh  =  {vh  e  Xh  :  (V  •  Vh,  qh)  =  0  Vqh.  G  Qh}  ■ 

Taylor-Hoocl  elements  (see  e.g.  [BS94,  G89])  are  one  common  example  of  such  a  choice  for  ( Xh,Qh ),  and 
are  also  the  elements  we  use  in  our  numerical  experiments.  Define  the  usual  explicitly  skew  symmetrized 
trilinear  form 

b*(u,  v,  w )  :=  -  (u  ■  Vi>,  w )  —  -(u  ■  Vw,  v). 
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2.2.  Discrete  differential  filters.  In  the  periodic  case,  for  <f>  £  L2(fi)  and  <5  >  0  fixed,  denote  the 
result  of  filtering  operation  on  </>  by  <j>,  where  4>  is  the  unique  solution  (in  X)  of 

— 52A(j)  +  <f>  =  (f>.  (2-6) 

We  let  G  :  L2(Q)  — >  X  denote  the  filtering  operation,  i.e.  <f>  :=  G(</>).  Following  Manica  and  Kaya- 
Merclan  [MKM06]  the  discrete  differential  filter  is  given  as  follows.  For  </>  £  L2(fl),  for  a  given  filtering  radius 
6  >  0,  Gh  :  L2(D)  — >  Xh,  (f)h  :=  Gh{f>)  where  <\>h  £  Xh  is  the  unique  solution  of 

<52(V4,V0  +  (4,0  =  ((/>,  vh)  Vvh£Xh.  (2.7) 

In  the  case  of  internal  flows  under  no  slip  boundary  conditions  we  must  use  the  discrete  Stokes  differential 
filter  to  preserve  discrete  incompressibility.  For  cj>  £  X* ,  6  >  0  given,  Gh  ■  X*  — >  Xh,  f>h  =  Gh{4>)  where 
{4>h,p)  G  Xh  x  Qh  is  the  unique  solution  of 

+  ((j>h,vh)  -  (p,  V-  vh)  =  (<j>,vh)  \/vh£Xh,  (2.8) 

(V-^,5)=0  Vq£Qh.  (2.9) 

The  Stokes  filter  results  in  an  exactly  divergence  free  filtered  velocity;  the  discrete  Stokes  filtered  velocity  is 
discretely  divergence  free. 

We  begin  by  recalling  from  [BIL04,  MKM06]  some  basic  facts  about  discrete  differential  filters. 

Lemma  2.1.  For  <j)  £  X,  we  have  the  following  hounds  for  the  discrete  filter 

fell  <  Ill'll ,  IIV0J  <  IIWH  and  ||V  x  4||  <  ||V4  . 

For  4>  £  X  and  Acj)  £  L2(fl) 

<52||V(^)-4)||2  +  H^-4112  <  C  inf  (62\\V(<t>-vh)\\2  +  \\<t>-vh\\2)+C64\\A(f>\\2. 

We  shall  assume  that  the  solution  to  the  NSE  that  is  approximated  is  a  strong  solution  and  in  particular 
satisfies 

u  £  L2{ 0,  T;  X)  n  L°°(0,  T;  L2(fl))  D  i4(0,  T;  X),  (2.10) 

p  £  L2(0,T;  Q),  ut  £  L2(0,T;  X*) ,  (2.11) 

and 

(ut,  v)  +  (u  ■  Vu,  v)  —  (p,  V  •  v)  +  v(Xu,  Vu)  =  (/,  v)  Vi>  £  X,  (2.12) 

(V  ■  u,  q)  =  0  Vq  £  Q.  (2.13) 

For  notational  clarity  let  v(tn+1^2)  =  v((tn+1  +tn)/2)  for  the  continuous  variable  and  vn+1^2  =  (vn+1  + 
vn)/2  for  both,  continuous  and  discrete  variables. 

2.3.  Deconvolution.  There  are  many  known  procedures  for  deconvolution,  e.g.,  [BB98].  The  minimal 
conditions  we  assume  throughout  are  that  the  (discrete)  filter  and  (discrete)  deconvolution  used  satisfy  the 
consistency  conditions  of  Stanculescu  [S07] . 

Assumption  Al:  Dh,  and  Gh  are  symmetric,  positive  definite  ( SPD )  operators  and  preserve  discrete 
incompressibility, 

Assumption  A2:  | \DhGh\ \c(L*->L2)  <  1  and  \\I  -  DhGh\\c(L^L^)  <  1, 

Assumption  A3:  (/  —  DhGh )  and  DhGh  are  SPD. 

These  have  been  proven  to  hold  for  van  Cittert  deconvolution  (next)  in  Stanculescu  [S07] .  Our  error  analysis, 
while  particular  to  a  differential  filter  and  the  van  Cittert  deconvolution  operators,  can  readily  be  extended 
to  more  general  ones  satisfying  the  above.  The  Nth  van  Cittert  deconvolution  operator  is  computed  by 
repeated  filtering.  It  can  be  compactly  given  as  follows.  (We  will  suppress  the  dependence  of  D  and  Dh  on 
the  parameter  N.) 

Definition  2.2.  The  Nth  order  van  Cittert  continuous  and  discrete  deconvolution  operators  as  D  and 
Dh  are,  respectively, 

N  N 

D<j>  :=  £(/  -  G)nct>  ,  and  Dh<j>  :=  ^(/  -  GhT4>  ■  (2-14) 

n=0  n=0 
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3.  Evolve  then  Filter  for  the  Navier-Stokes  equations.  As  a  first  step,  we  consider  discretization 
by  the  FEM  in  space,  Crank-Nicolson  (CN)  method  in  time  with  an  added  filtering  step  and  no  deconvolution 
or  relaxation  (i.e. ,  Step  1  and  Step  2a  of  Algorithm  1.1).  If  no  slip  boundary  conditions  are  imposed,  let 
the  filter  be  the  discrete  Stokes  filter  (2. 8), (2. 9),  and  if  periodic  boundary  conditions,  the  simple  discrete 
differential  filter  (2.7),  denoted  by  G  and  Gh- 

Algorithm  3.1.  [Evolve  then  Filter  for  NSE] 

Step  1:  Given  ulf  find  £  Xh,  Ph+1  £  Qh  satisfying 


,w 


n+1 


-  Uh 


At 


,vh )  +  h*( 


w 


n+1 


+  Uu  w 


n+1 


-,Vh)  +  +(V 


W 


n+1 


-,V+)-(K+1,V 


Vh) 


=  (fn+1/2,vh),  for  all  vh  £  Xh, 
(V  •  w%+1,qh)  =  0,  for  all  qh  £  Qh. 


(3.1) 


Step  2:  Filter  z++1  to  give  v% 


n+l  satisfying 
S2(Xu 


l+\Vvh)  +  {unh+\vh)  =  {w 


imposing,  in  the  non-periodic  case,  (V  •  u 


n+1 
h 
n+1 
h  ’ 


n+1 

h 


,  Vh)  Vvh  £  xh  . 
qh)  —  0,  ^  qh  £  Qh  • 


The  temporal  consistency  error  and  associated  forecasted  global  error  in  the  basic  algorithm  (3.1)  is: 


Temporal  Consistency  Error  =  0(At3  +  52),  Global  Error  =  0(At2  + 


S2 

At 


+  Spacial  Error) , 


We  turn  to  stability. 

Lemma  3.2.  Let  u =  vihh-  We  have 

2^2||V<||2+  ||<||2+  |K-<||2  =  IK||2. 

Proof.  Consider  the  discrete  Stokes  filter  (the  proof  is  the  same  in  the  periodic  case).  Recall  that 

<S2(VuJ*,  Xvh)  +  (u%,vh)  =  «,vfc)  Vvh  £  Vh. 


Set  Vh  =  uf-  This  gives 

^||v<||2  +  |K||2  =  «,<). 

The  polarization  identity  (w%,  u%)  =  5||i+||2  +  g||«/(|i2  —  ^\\uh~wh\\2  gives,  after  rearrangement,  the  result: 

IKII2  =  2<52||V<||2  +  IKH2  +  |K  -  <||2. 


Next  we  prove  a  strong  energy  equality  and  associated  strong,  unconditional  stability  property. 
Proposition  3.3.  The  approximate  velocity  u][+1  given  by  the  Algorithm  3.1  satisfies  the  energy  equality 


\  I K+1 1 12  +  1 1  V<+1 1 12  +  \K+1  ~  <+1  II2' 

n= 0  ' 


/!|v( 


wh  +  1  +  <  u,2 


=  2  U 


,  0 1 1 2 


n+1/2  wh 


n+1 


Wh  -Yu 


n— 0 


and  the  stability  bound 
1 


n— 0 


cn2+A<i:(^iiwriii2+;V|i«r-<+'ii2+^iiv(“''- 


—  | 
2AC 


<  xllA 


0  1 1 2 


r 


+1/2  1 1 2 


n— 0 
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(3.2) 


Proof.  In  Step  1  set  Vh  =  (w£+A  +  v%)/2.  This  gives 

^(IK+1II2-IKII2)  +  H|V(^T±U" 


ft.'ill2  _  (fn+ 1/2 


n+1  „,n 


We  use  the  stability  equality  of  Lemma  3.2 


n+ 1 1 1 2 


1 1 wh  '  IT  —  2i52||Vw)(+1||2  -f:m 
in  the  first  term  in  the  LHS  of  (3.2).  Rearranging  gives 


K+1ll2  +  IK+1-<+1ll2 


2A  ^ 


.n+1  1 12 


1 


-  IKH2)  +  (  ^  KK+r  +  ^IK+1  -  w 


n+ 1  n+l||2 


'l|V( 


<+  +  <\||2 


_  (yn+1/2 


n+1  i  „.n 


+ 


Summing  this  establishes  the  energy  equality.  Using  the  Cauchy-Schwarz- Young  inequality  on  the  RHS, 
subsuming  one  term  into  the  LHS  gives 


2Af( 


,n+ 1 1 1 2 


(52 


1 


-  IK II  )  +  ^IIVK  II  +  2^IIK+i  -  <11  +  2  K( 

<^iir+1/2n2 


I/M  +  u 


n+ 1  |  n 


h  \||2 


Summing  over  the  index  n,  the  global  stability  estimate  follows. 

□ 

The  method  is  thus  stable.  The  viscous  and  numerical  dissipation  in  the  method  are  respectively 


Viscous  dissipation  :=  i/||V( 


Wh+  +K\|i2 


Numerical  dissipation  :=  —  ||Vi+|2  +  — — 


n+1  n+li|2 


A t"  "  “7l"  '  2 A t"“h 


Uu  -w 


The  first  term  in  the  numerical  dissipation  is  the  same  as  induced  by  the  explicit  method  plus  simple  filtering. 
The  second  term  resembles  the  numerical  dissipation  in  the  backward  Euler  method. 

4.  Evolve,  Filter,  Deconvolve  and  Relax.  Consider  discretization  of  the  Navier-Stokes  equations 
by  the  FEM  in  space,  CN  method  in  time  with  an  added  filtering  step  including  deconvolution  and  time 
relaxation.  To  give  concrete  estimates  specific  choices  of  deconvolution  operators  and  filters  must  be  made. 
In  this  section  we  choose  the  (discrete)  van  Cittert  deconvolution  operator  and  (discrete)  Stokes  filter. 

Algorithm  4.1.  [Evolve,  Filter,  Deconvolve  then  Relax  for  NSE]  Let  the  filter  he  the  discrete  Stokes 
filter  (2.8),  (2.9).  Choose  x  with  0  <  x  <  1. 

Step  1:  Given  u find  <+1  G  A+  K+1  ^  Qh  satisfying 


wn+ 1 _ ,,n 

(^AK  +  K 


^+1 


w^-+uh  wh 


,71+ 1 


n+1 


’  2  ’  '  v  2 
=  ( fn+1/2,vh ),  for  all  vh  G  Xh, 


yvh)-{pnh+1y  -vh) 


(V 


n+1 


qh)  =  0,  for  all  qh  G  Qh. 


Step  2:  Compute  Dh  »  then  relax: 


(4.1) 

(4.2) 


uh+1  =  (1  -  xK+1  +  xDh  (<+1h) 


We  next  analyze  the  methods  enhanced  stability. 


(4.3) 

(4.4) 


Lemma  4.2.  Let  Uh  =  (1  —  x)wh  +  X^h  ( u>hh )•  Then 

I  Ml2  -  IKII2  =  X(2  -  X)((I  ~  DhGh)wh,wh)  +  x2((/-  DhGh)wh,DhGhwh), 

IKII2  -  IKII2  =  -\\uh  -  wh\\2  +  2x({I  -  DhGh)wh,wh) . 

If  Assumptions  Al,  A2  and  AS  hold  then 

IKII  <  IKII-  (4.5) 

Proof.  The  first  equality  is  proven  in  Theorem  1.2.  For  the  second,  take  the  inner  product  of  = 
(1  -  xK  +  XDhGhwh  with  wh.  This  gives 

( uh,wh )  =  (1  -  x)(wh,wh)  +  x{DhGhwh,wh)  =  |K|| 2  -  x{{!  ~  DhGh)wh,wh). 

Applying  the  polarization  identity  to  the  LHS  we  obtain 

^ IKH 2  +  t,IKII2  -  *IK  -  Kl2  =  ( uh,wh )  =  IKH2  -  x({i  -  DhGh)wh,wh) , 

from  which  (4.4)  follows.  For  (4.5),  note  that  since  ||KG/,||  <  1  and  0  <  x  <  1, 

IKII  <  {t-x)\\wh\\  +  x\\DhGhWh\\  <  ((1  -X)  +x\\DhGh\\)  IKII  <  |KII- 


□ 

Next  we  prove  an  energy  equality,  unconditional  stability  and  give  the  precise  formula  for  the  numerical 
dissipation  in  the  algorithm. 

Proposition  4.3.  [Stability  with  Deconvolution  and  Relaxation] Suppose  Assumptions  Al,  A2  and  AS 
above  hold.  Algorithm  f.l  satisfies  the  energy  equality 


1, 


Z+l  1 1 2 


-pAf^||V( 


W 


n+1  I  „.n 


*h \||2 


n— 0 


n— 0 


-  DhGh)wnh+\  <+1)  +  ^((/  -  DhGh K+1,  DhGhwnh+1) 


1 

=  vll  u 


0  1 1 2 


,„”+l  4-  nn 

+  A tJ2(.fn+1/2,  h  2  K 


n— 0 


and  t/ie  stability  bound 


4+1ll2 


pAt^||V( 


w 


n+1 


xMl|2 


n= 0 


+At^ 


n— 0 


^K^((/  -  KKK+\  <+1)  +  ^((/  -  KKK+1,  KK<+1) 
<IKII2  +  —  £|ir+1/2||2. 


z/ 


n=0 


^(IK+1II2 


!)  +  Hlv( 


n+l  |  n 


Proof.  In  Step  1  in  Algorithm  4.1  set  Vh  =  (i++1  +  t/^) /2.  This  gives 

1 


'■'"2  1  +[IK+iII2-IK+iII2]  =  (/"+1/2.!+ 


n+1  I  „,n 


-)•  (4-6) 


We  use  Lemma  4.2  in  the  bracketed  term  in  the  LHS  of  (4.6).  Rearranging  gives 


^(IKKI2 


!)  +  H|v( 


w 


n+1  I  „,ra 


Ami2 


^^((J  -  DhGh)wp\ <+1)  +  ^((/  -  KKK+1,  DhGhwnh+1) 


_  (j?n+ 1/2  Kl 


n+1  I  „.n 


-)• 
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Summing  this  proves  the  energy  equality.  Using  the  Cauchy-Schwarz- Young  inequality  on  the  RHS,  sub¬ 
suming  one  term  into  the  LHS  and  summing  over  the  index  n,  proves  the  global  stability  estimate. 

□ 

The  numerical  dissipation  in  Algorithm  4.1  is  exactly  described  by  the  bracketed  term  in  the  energy 
equality: 


Numerical  dissipation  := 


X(2  ~  x) 

2A  t 


((/  -  DhGh)w\ 


n+1  n+1 


)+2Kt{{I~  DhGh)wnh+1,DhGhwnh+1) 


5.  Error  Analysis  of  the  Time  Relaxation  Approximation.  In  this  section  we  present  a  detailed 
error  analysis  for  the  approximation  scheme,  Algorithm  4.1  incorporating  both  (differential)  filtering  and 
(van  Cittert)  deconvolution.  We  focus  our  attention  on  the  simplest  differential  filter,  G  defined  via  (2.6) 
and  its  discrete  counterpart  Gh  defined  by  (2.7).  Lemma  2.1  presents  the  error  estimate  for  the  error  between 
t p  and  DhGh{<t> )  :=  Dofih  for  the  0th  order  deconvolution  operator.  For  van  Cittert  deconvolution  operators 
(Definition  2.2)  we  have  the  following  result. 


Lemma  5.1.  [LMNR.06]  For  smooth  <f>  the  discrete  Nth  order  approximate  deconvolution  operator  satis¬ 
fies  for  0  <  s  <  N 


H  -  DhGh<t> II  <  C\  52s+2  ||0||ij2»+2 


C2  ( 5hk  +  hk+1 ) 


'JV+1 

E  \°n^)\ 


i 


(5.1) 


The  dependence  of  the  +"(/>)  |fc+i  terms  in  (5.1)  upon  the  filter  radius  <5,  for  a  general  smooth  function 
</>,  is  not  fully  understood.  In  the  case  of  <j>  periodic  the  \Gn((f)\k+i  are  independent  of  5.  Also,  for 
satisfying  homogeneous  Dirichlet  boundary  conditions,  with  the  additional  property  that  AJ</>  =  0  on  dfl 
for  0  <  j  <  [^+]  —  1)  the  \Gn(<j))\k+i  are  independent  of  6,  see  [LMNR06],  [L07]. 

As  mentioned  above,  Taylor-Hood  approximating  elements  are  a  common  choice  for  ( Xh,Qh )  and  cor¬ 
respond  to  k  =  2  in  (5.1).  For  these  approximating  elements  we  have  from  Lemma  5.1  the  following  corollary. 


Corollary  5.2.  Suppose  <j>  £  Hq(Q)  n  if4(fl).  Suppose  the  order  of  deconvolution  is  N  =  1  and 
( Xh,Qh )  are  chosen  to  be  the  Taylor-Hood  elements.  We  have 

U  -  DhGh(t>h II  <  Ci  <53||<(>||3  +  C2  ( 6h 2  +  h3)  ||0||3  .  (5.2) 

Proof.  This  follows  from  the  previous  lemma  by  taking  s  =  1/2,  k  =  2,  N  =  1  and  thus  —1  =  0. 

We  have  then  ||<^||3  <  C||^||3  with  constant  independent  of  5  and  h. 

□ 

Extending  Corollary  5.2,  we  will  make  the  following  assumption. 

Assumption  DG1:  The  \Gn((j))\k+i  terms  in  (5.1)  are  independent  of  <5  and 

\\fi  -  DhGh<f> II  <  CiS2N+2  U\\H2N+2  +  C2  (6hk  +  hk+1)  HWk+i.  (5.3) 


With  tn  =  nAt,  n  =  0,1,2 ,...,NT,  T  :=  NTAt,  and  dtfn  :=  ( f(tn )  —  f(tn  l)/A t,  we  introduce  the 
following  discrete  norms: 


v  oo.fe  :=  max  u  fc. 

0<n<NT 

/  Nt 


1/m 


IIK/2|||oo,fe  :=  max  ||vn  1/2||fc, 

1  <n<NT 

f  Nt  ^ 

II  1^1/2 1  ||m,fc  :=  Ell^_1/2H“At 


1/m 


\n— 0 


We  denote  time  level  averages  by 


w  := 


wn  +  un  1 


and  wn  := 
10 


wn  +  wn  1 


(5.4) 


To  begin  the  analysis  we  rewrite  Algorithm  4.1  in  the  following  form. 

Algorithm  4.1  restated.  Assume  that  the  filtering,  deconvolution,  and  relaxation  parameters  S,  N, 
and  x  are  given.  Then,  for  n  =  1,2, ,  Nt,  find  w%,  DhGh{w %),  u ^  G  Xh,  p7/  G  Qh,  such  that 

O h,vh)  +  A tb*(w%,iv%,vh)  -  At  {pi,  V  •  vh)  +  Atv(Vw%,  Vvh) 

=  (up\vh)  +  At  (fn~1/2,  vh) ,  Mvh  G  Xh  ,  (5.5) 

(X-wZ,qh)  =  0,Vqh£Qh,  (5.6) 

<  =  (1  -xX  +  xDhGh{wl).  (5.7) 

As  the  spaces  Xh  and  Qh  satisfy  the  usual  inf-sup  condition,  we  can  equivalent  consider  the  problem: 
For  n  =  1,2, ... ,  Nt  find  w7/,  w%h ,  u7/  G  Vh  such  that 

(wh,vh)  +  A  tb*(ivh,Wh,vh)  +  Atv{S7wh,Vvh)  =  (v%fi%vh)  +  At  (/”~1/2,  vh),  Vvh  G  Vh, 

S2(V  ■  w%h,  V-vh)  +  (w%h,  vh)  =  (w%,  vh)  ,  Mvh  G  Vh  , 

<  =  (1  -  X)wh  +  xDhKh  ■ 

We  first  establish  computability  of  the  approximation  at  each  time  step. 

Lemma  5.3.  For  the  approximation  scheme  (5.8)-(5.10)  we  have  that  w%,  whhi  uh>  exist  ad  each  time 
step. 

Proof.  The  existence  of  a  solution  w7/  to  (5.8)  follows  from  the  Leray-Schauder  Principle  [Z95] .  Specifi¬ 
cally,  with  A:Vh~>  Vh,  defined  by  y  =  A(w) 

( y,v )  :=  —Atb*((w  +  up1)/ 2,  (w  +  up1)/ 2,v)  —  Atv(S7(w  +  up1)/ 2,  Vv) 

+  (up\v)  +  At(/n_1/2,  v)  , 

the  operator  A  is  compact  and  any  solution  of  w  =  s  A(w) ,  for  0  <  s  <  1  ,  satisfied  the  bound  ||w;||  <  7, 
where  7  is  independent  of  s.  Existence  and  uniqueness  of  w7/,  follows  directly  from  the  Riesz  representation 
theorem.  Existence  and  uniqueness  of  u ^  follows  for  w7/  and  w7/,  from  the  definition  of  Dh- 
□ 

To  establish  the  optimal  asymptotic  error  estimates  for  the  approximation  we  need  to  assume  that  the 
true  solution  is  more  regular  than  that  given  by  (2. 10), (2. 11).  Specifically,  assuming  2 N  +  2  >  k  +  1, 

uG  Loo(0,T;W4fc+1(9))nLoo(0,T;iJ2Ar+2(9))n  (5.11) 

7^(0,  T;  Hk+1(n))nH3(0,T-,  L2{Q))  n  W42(0, T;  H1^)), 

pG  L°°{0,T-HS+1(CI)),  and  /  G  H2(0,T- L2(Q)) .  (5.12) 


(5.8) 

(5.9) 
(5.10) 


For  the  errors  en  :=  u(tn)  —  u7/,  and  en  :=  u(tn)  —  w7/,  we  have  the  following. 

Theorem  5.4.  Foru,  p,  and  f  as  described  by  (5. 11), (5. 12),  satisfying  (2. 12), (2. 13),  andu7/,  w7/  given 
by  (5.5)-(5.7)  we  have  that  for  At  sufficiently  small 


\u  Uh 1 1 1 00 ,0  +  |||w  -  w/t|||oo,o  <  F(At,  h,  6,  x)  +  Chk+1\\ |u|||oo,fc+i , 

(^AtJ2\\V(un~1/2  -  {wr/l+up1)/2)\\2^  <  F(At,h,5,x)  +  CP/2(At)2\\ Vutt||2,0 

+  CP/2hk\\\u\pk+1  ,  for  1<  l  <Nt. 


(5.13) 

(5.14) 
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where 


F{At,h,6,x)  ■■=  Cv~1/2  (hk+1/2\\\u\\\lk+1  +  *fc+1/2|||Vu|||!i0  +  /is+1|||pi/2|||2,s+i) 

+  CV1/2  hk\\\u\\\2,k+i  +  Chk+1\\ut\\2,k+1+ Cv-1hk\\\u\\\00,k+1 
+  Cx  (w3/2|||Vu|||^)0  +  (At)-1/2)  (Aty1/2  hk+1  \\\u\\\2,k+i 

+  Cx  (v  3^2  ||  I V  II  ^o>0  +  (At)  (At)  !/2  (<52W+“  +  (5/lfc  +  /lfc+1)  (|||u|||2,2iV+2  +  HM||2,fc+l) 
+  C  (At)2  ^||utit||2,0  +  ||/tt||2,0  +  "1/2  |Vutt||2)0 

+  ^  1/2||Vwtt||4,o  +  v  1/2|||Vw||||i0  +  v  1/2|||VUl/2||||,o)  • 

Proof.  At  time  t  =  (n  —  1/2) At,  u  given  by  (2.12)-(2.13)  satisfies 

{un-un~\vh)  +  ^  v(A7un  ,  Vvh)  +  Atb*(un  ,un  ,  vh)  -  At  {pn~1/2,  V  •  vh) 

=  A t(fn,vh)  +  A tlntp(un;vh)  ,  (5.15) 

for  all  Vh  G  Vh,  where  Intp(un\Vh),  representing  the  interpolating  error,  denotes 

Intp(un;vh)  =  ^( un  —  un_1)/At  —  u"-1^2 ,Vhj  +  v(yun  —  Vwn~1/2  ,  Vvh) 

+  b*{un,un,vh)  -  ^(n""1/2,!*""1/2,^) 

+  (fn-1/2-fn,vh).  (5.16) 

Subtracting  (5.8)  from  (5.15),  we  have  for  en  =  un  -  w%,  en  =  un  -  ujj, 

{en-en~\vh)  +  At  v(y  (en  +  e"_1)/2 ,  Vuj,)  =  —  At  b*{un  ,  un  ,  vh)  +  Atb* (wnh,wnh,vh) 

+  A t(pn_1,/2,V  •  Uft)  +  At  Intp(un-,Vh)  ,  (5-17) 

for  all  vh  G  14.  Let  E/n  G  Vh,  en  =  un  -  w%  =  (un  -  Un)  +  (Un  -  w%)  :=  An  +  Fn,  and  en  = 

un  -  u%  =  ( un  -  Un)  +  {Un  -  ul)  :=  A"  +  En. 

Similar  to  the  notation  defined  in  (5.4),  we  use  Fn  :=  ( Fn  +  En~1)/2.  With  the  choice  Vh  =  Fn,  and 

using  (V  •  Fn,qh)  =  0,  \/qh  G  Qh,  equation  (5.17)  becomes 

i(||FnH2  -  ||£,r!-1||2)  +  At ^|| VAm||2  =  -(A”  —  An~l  ,Fn)  -  At^(VA”,VF”) 

-  A tb*{un  ,un  ,  Fn)  +  A tb*{wl,wl,Fn) 

+  At(p”-1/2  -qh,V-  Fn )  +  At  Intp(un ;  Fn)  .  (5.18) 

Next  we  estimate  the  terms  on  the  RHS  of  (5.18). 


1  An  _  A”-1  1 

(An  —  An~1,Fn)  <  —  At || - — - ||2  +  -At||F"||2 


1 


1 


=  2At L [mLiAtdt  dn  +2At||F 


in||2 


/*"- 1 
rtn 


i 


f  ||At||2dt  +lAt(||^||2  +  ||i?"-1||2)  .  (5.19) 

Jtn~ 1  Z 
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(5.20) 


j/(VA”, VF”)  <  —  ||VF"||2  +  C'i/||VAn||2. 

We  rewrite  b*(un  ,un  ,  F”)  —  b*(w^,u^,  F”)  as 

b*  (un  ,  un  ,  F”)  -  b*(w^,  iv^,  F")  =  b*{un  ,  un  ,  F”)  -  b*{w%,un,Fn) 

+  b*(Wfr,un,Fn)  -  b*{wnh,wnh,Fn) 

=  b*((en  +  e”_1)/2  ,un  ,  F")  +  b*(w^,  (sn  +  en~1)/2,  F”) 

=  b*(kn  +  Fn  ,un  ,  F”)  +  b*(Wfr,  A”  +  F”  ,  F”) 

=  b*(An,un,  Fn)  +  b*(Fn  ,  un  ,  F”) 

+  b*(w%,  A"  ,  F”)  +  ft*(i»J*,  F”  ,  F”) .  (5.21) 

Using  b*(u,v,w)  <  C(U)-\/||w||  ||Vw||  ||Vw||  ||Vw||,  for  u,  v,  w  €  X ,  and  Young’s  inequality,  we  bound 
the  terms  on  the  RHS  of  (5.21)  as  follows. 


b*(An  ,un  ,  Fn)  <  Cy||A"||  || VAn||  ||Vw"||  1 1 V Fn \ 

<^I|VF”||2  +  C  v~l  ||A"||  ||  VA”  ||  ||Vw”||2 


(5.22) 


b*(Fn  , un  ,  F”)  <  C||F”||1/2||VF”||3/2||Vu”|| 

<^||VF”||2  +  C u~3  ||Vu”||4  ||F"||2 

<  ^||VFn||2  +  C is~3  ||Vu”||4  (||F”||2  +  ||Fn_1||2)  (5.23) 

b*(w%,An,  F”)  <  C||V<||  ||VA"|  ||VF”|| 

<  ^  l|VFn||2  +  Cv -1  ||V^||2  || VA” || 2  (5.24) 

6*(<,F”,F”)  =  0  (5.25) 

(pn~1/2  ~  qh,  V  •  F”)  <  | Ip""1/2  -  ®,||  ||  V  •  FI 

<  ^l|VF”||2  +  Cv-1  lb”"1/2  -  9b|2  •  (5-26) 

With  the  bounds  (5.19)-(5.26),  (5.18)becomes 

i(||F”H2  -  ||F”-1||2)  +  A^||VF”||2 

<  CAt(l  +  v~3\\Vunf)  (\\Fn\\2  +  IIF”-1!!2)  +  CYAt||  VA"||2 
+  CV-1Af||Vw£||2||VA”||2  +  CV_1Af  ||  Vm"||2||A"|  ||VA”|| 

tn 

+  C  f  ||At||2  dt  +  Cv~1At  |b"-1^2  -  Qh\\2  +  At  \Intp(un;  F”)| .  (5.27) 

Jtn~x 

As  u ^  and  w ^  are  connected  through  the  filter- deconvolve-relax  equation,  we  next  use  that  equation  to 
obtain  a  relationship  between  ||F”||  and  ||F™||.  The  true  solution  u(-,tn)  =  un  satisfies 

«"  =  (!-  xK  +  xDhGhun  +  X(I  ~  DhGh)un  .  (5.28) 
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Subtracting  (5.7)  from  (5.28)  yields 


e™  =  (1  -  x)en  +  xDhGhen  +  X(J  -  DhGh)un  (5.29) 

i.e.,  En=  (1  -  X)Fn  +  xDhGhFn  -  X(J  -  DhGh)An  +  X(J  -  D,,Gh)un 
>  \\En\\  <  (1  -  x)\\Fn\\  +  x\\DhGh\\  ||F”||  +  X||(J  -  DhGh) A"||  +  X||(J  -  DhGh)un\\ 

\\En ||  <  ||F"||  +  x\\(I  ~  DhGh)An\\  +  x\\(I  -  DhGh)un\\  (5.30) 


Squaring  both  sides  of  (5.30)  and  simplifying  we  obtain 


1 1  En  1 1  2  <  ||Fl2  +  Af||Fl2  +  2x2(l  +  (At)'1)  \\(I—  DhGh)An\\2  +  2X2  (l  +  (At)-1)  \\(I-DhGh)unf  . 

(5.31) 

Substituting  (5.31)  into  (5.27)  and  summing  from  n  =  1  to  l,  (assuming  that  ||F°||  =  0,  i.e.  u°  £  Xh) 
we  obtain 


yin  ||  2 


7,n\\2 


||F‘||2  +  z,Af^||VF 

n=l 
l 

—  C  (1  +  1,  3  l|Vttn||4)  H-F 

n— 0 

Z  l 

+  CvAt,J2  ||VA"||2+  C ^-1  At  ^2  II Vu"||2  ||A”||  || VA 

n= 1  n= 1 

+  C'^1At^||V<||2||VA”||2  Ecj^f  llA*ll2^ 

n=l  n— 1  ^n  — 1 

1-1 

+  C  At  ^2  (l  +  II  Vwn||4  +  (At)-1)  (1  +  (A t)-1)  x2\\(I  ~  DhGh)A 

n— 0 
l-l 

+  C  At  ^  (l  +  v~3  || VMn||4  +  (At)-1)  (1  +  (At)-1)  X2||(I  —  DhGh) 

n— 0 

+  At  ^2  C  ^-1||p”-1/2  —  qh\\2 

n—  1 
l 

+  2 At  J2  \Intp(un;  Fn)  \ . 


n  II  2 


n= 1 


The  terms  on  the  RHS  of  (5.32)  can  be  further  simplified  as  follows. 


(5.32) 


Cv  At  ^2  ||  VAn||2  <CuAt^\\XAnr  <  Cv  At  <Cvh™  |||«|||^fc+1 


-  \  llT7An||2 


L.2k\„.n\2 


n— 0 


n=0 
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(5.33) 


For  the  next  term 


CV_1  At  ^2  l|Vwn||2  ||A"||  ||VAn|| 

n= 1 

Z 

<  Cv^At  J2  (l|A”||  ||VA"||  +  HA”-1!!  IIVA"-1!!  +  HA"-1)!  II VAn||  +  ||A"||  UVA"-1!!)  || Vu" 

n— 1 

/  l  l 

<  Cv-1  h2k+1  At  K\l+1  II Vm"||2  +  At.^2  |w"|fc+1|M"_1|fe+1 1| Vm"||2 


+  At^2\ 


<  Cv~l  h2k+1  A t.J2 


Z^\  I&+1  _r  , 

n— 0  n= 0 


,n- 1 1 2  ||v7-ti||2 


A ||Vu" 


<Ci/-I/l2fc+1(||H||tfc+1  +  |||v«|||to)  • 

Using  the  boundedness  of  v  At  J2n=i  l|Vt&£||  (Proposition  4.3) 
i  i 

Cv-1  At£  ||V«;£||2||VA12  <  Cu^AtJ^  l|V^||2-/i2fc(|K||2fc+i  +  IK^IlLi) 

n— 1  n— 1 

<  Cv~2h2k  HMII^fc+i  • 


Cj2fl  \\Atfdt  <  cj^  f  h2k+2\\ut\\2dt  <  Ch2k+2\\ut\\lk 
Jt "-1  “i  .it"-1 


AtY^Cv-^-^-qnW2  <  Cv-'At^h^Wp^'Ys+i  <  Cv~x  h2s+2  |||p1/2||||,s+1 .  (5.37) 

n=l  n—1 

From  Lemma  5.1  and  Assumption  A2,  and  assuming  At  <  1 
l- 1 

<7  At  5^(1  +  ^-3  || Vm”||4  +  (At)-1)  (1  +  (At)-1)  X2||(J  -  DhGh) A"||2 

71=0 

Z-l 

<  (7  At  ^  (v~3  ||Vu"||4  +  (At)”1)  (At)-1  X2\\I  -  DhGh\\2  \\Anf 

71=0 

Z-l 

<  CAtYl  (v~3  ||Vu"||4  +  (At)-1)  (At)-lx2h2k+2\\un\\l+1 

71=0 

<  C'X2  (^-3|||VU|||l,o  +  (At)-1)  (At)-1  h2k+2  HMIll.fc+i  •  (5-38) 


GAt^2  (l  +  ^-3||Vi7n||4  +  (At)-1)  (1  +  (At)"1)  x2\\(I  ~  DhGhju" 


<  CAt^  {V~Z  l|Vwn||4  +  (At)-1)  (At)-1  (<54JV+4  +  S2  h2k  +  h2k+2)  x2  (\\un\\2N+2  +  \\un\\2k+1) 

71=0 

<  C x2  (y  3  II  IVulH^o  +  (At)  x)  (At)  1  (S4N+i  +  52  h2k  +  h2k+2)  (\\\u\\\2  2n+2  +  HMIli.k+i)-  (5-39) 
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As  given  in  [ELN07]  the  interpolation  error  term  2A t  Y^n= 1  \Intp(un;  Fn)\  in  (5.32)  can  be  bounded  as 
i  i  i 

2A  t  J2  \Intp{un\  Fn)\  <  A  tC^  ( \\Fnf  +  ||^n_1||2)  +  (ex  +  e2  +  e3)At  v  ^  || VFn\\2 

n= 1  n= 1  n= 1 

+  C  (At)4  (WutuWtfi  +  II /ct II 2,0  +  ^l|Vw«||i)0 

+  ^“1||Vutt||t0  +  ^IIIVwHlto  +  ^-'llIVtii/alHto)  •  (5-40) 

Combining  (5.33)-(5.40),  equation  (5.32)  simplifies  to 

i 

\\Fl\\2  +  At^v\\S7Fn\\2 

n— 1 

l 

<  CAt  J2  (!  +  ^_3|| Vm"||4)  ||F"||2 

n= 0 

+  Cv  1  (/l2fe+1|||u||||fc+1  +  ^2fe+1|||VM||||>0  +  fr“S+2|||Pl/2|||2,s+l)  +  C'zy^2fe|||w|||2!fc+1 
+  C  h2k+~  \\ut\\\ik+i+  Cv  2 /i2fe|||«|||^0)fe+i 
+  C\2  (v~3  |||Vu|||^)0  +  (At)^1)  (Aty1  h2k+2  \\\u\\\lk+1 

+  C x2  (v  3  || | Vi*|  11^0,0  +  (At)  4)  (At)  1  (£4Ar+4  +  S2  h2k  +  h2k+2)  (\\\u\\\lt2N+2  +  II  Mlll.fc+i) 

+  C  (At)4  (||wttt||2,o  +  ll/**lll,o  +  HIVuttlli.o 

+  ^1||VMtt||40  +  iz-'IIIVulllto  +  z/-1|||Vu1/2|||4)o)  •  (5.41) 

Hence,  with  At  sufficiently  small,  i.e.  At  <  C(1  +  ^_3|| | Vit| H^,  0)_1,  from  the  discrete  Gronwall’s 
Lemma  [HeRa],  we  have 

i 

||F'||2  + At  ^V||VF"||2 

n= 1 

<  Cv-1  (h2k+l\\\u\\\4Ak+1  +  /i2fc+1|||Vu|||4i0  +  /i2s+2||bi/2llli,s+i)  +  c,^/i2fc|||w|||2fc+i 
+  C h2k+2  \\ut\\%k+i  +  Cv  2 /i2fe|||u|||^o  fc+i 

+  Cx 2  O'-3  |||Vm|||^>0  +  (At)^1)  (At)~4  h2k+2  \\\u\\\lk+1 

+  C x2  (v  3  IIIVuIH^q  +  (At)  4)  (At)  1  (c>4Ar+4  +  S2  h2k  +  h2k+2)  (|||w|||2 )2iv+2  +  II  1^1  II l,fc+i) 

+  C  (At)4  (\\uttt\\2,o  +  II/** lll,o  +  *^ll  V "**** II l,o 

+  ^_1||Vwtt||to  +  ^"'lllVulHto  +  I'-'lllVui/alllto)  ■  (5-42) 

The  estimate  given  in  (5.13)  for  |||u  —  Wh  | ||oo,o  then  follows  from  the  triangle  inequality  and  (5.42).  The 
estimate  for  |||u  —  *t/i|||oo,o  follows  from  the  estimate  (5.30),  the  triangle  inequality,  and  the  estimate  for 
|||w  —  Wh  1 1 1 oo,0 -  To  obtain  (5.14),  we  use  (5.42)  and 

||V  («n-1/2  -  «  +  <-1)/2)  II2  <  ||V(wn-1/2  -  un) ||2  +  || VA”||2  +  || VF"||2 

<  l|V««||a*  +  Ch2k\un\2k+1  +  Ch2k\u”-4\2+1 

+  ||Vi>||2. 

□ 

For  Taylor-Hood  approximating  elements,  i.e.  k  =  2,  s  =  1,  we  have  the  following  asymptotic  estimate. 
Corollary  5.5.  Under  the  assumptions  of  Theorem  5.4,  with  6  =  Ch,  At  =  Ch,  N  =  1,  and 
(Xh,Qh)  given  by  the  Taylor-Hood  approximation  elements,  we  have 

|||«  -  u>/,|||oo, o  +|||«  -  Ufcllloo.0  +  fi/At  ^IIVK*-1/2  -  K  +  <-1)/2)||2 
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<  C  ((At)2  +  h2)  . 

(5.43) 


The  filter-deconvolve-relax  step  is  more  than  a  simple  perturbation  of  the  usual  FEM  at  each  time  step. 
Computationally  its  positive  influence  on  the  stability  of  the  approximation  algorithm  is  very  apparent. 
Mathematically  to  see  this  increased  stability  we  need  to  consider  a  different  norm  than  the  L 2  (in  space) 
norm.  Define  the  (mesh  dependent)  weighted  norm  and  inner  product 

(' v,w)idg  ■=  ((I  -  DhGh)v,w) ,  and  \\v\\ idg  ■=  \/{v,v)IDG  ■  (5.44) 

From  Figure  1.1  we  see  that  this  norm  has  the  property  of  measuring  high  frequency  components  of  a  func¬ 
tion.  Estimate  (5.45)  shows  that  the  high  frequency  components  of  w^,  where  spurious  oscillations  would 
concentrate,  are  reduced  in  u%,  i.e.  |K||jdg  <  HKIKg-  Estimates  (5.46),  (5.47),  establishes  a  relation¬ 
ship  between  the  L 2  errors  of  u l  and  w^,  and  the  high  frequency  components  of  the  error  in  uj)  and  w%,  i.e. 
|| un  —  Ufl\\iDG  and  IK  —  w%\\idg,  respectively.  We  remark  that  the  terms  in  the  brackets  on  the  RHS  of 
(5.46)  and  (5.47)  represent  consistency  error  terms  for  the  filter- deconvolve  operation. 

Theorem  5.6.  Under  the  assumptions  of  Theorem  5.f,  for  n  =  1,2,...,  Nt,  0  <  <5  <  1, 

IKII  2idg  =  \\<\\]DG-x£-x)W-DhGh)wnh\\2 -X2  {{I-DhGh)wnh,  {I  -  DhGh)DhGhwl)  ,  (5.45) 

IK  -  <ll2  <  IK  -  Kll2  -  fx( i  -  x)IK  -  <11 2idg 

+  X2\\(I-DhGh)un\\2+^x(l  +  ^X)\\un\\jDG  +  x2((I-DhGh)un,DhGhun)  .  (5.46) 

IK  -  <11 2 dg  <  IK  -  wl\\]DG  -  X\{I-  DhGh)(un  -  <),(/  -  DhGh)DhGh(un  -  <)) 
-\x{l-xW-DhGh){un-wl)\\2 

+  [2X(1  +  x)IK  -  DhGh)un\\ 2  +  x2||(/  -  DhGh)un \\jDG]  .  (5.47) 

Proof.  We  note  that  from  Assumption  A3,  ((/  —  DhGh)v ,  DhGh  v)  >  0,  and 

((/  —  DhGhfv ,  (/  —  DhGh)DhGh  v)  >  0.  To  simplify  notation,  in  the  proof  we  use  e  :=  en  =  u(tn)  —  u'f 
and  s  :=  en  =  u{tn)  —  wjf  and  IDG  to  denote  I  —  DhGh •  From  (5.7),  taking  the  inner-product  of  both  sides 
with  respect  to  (I  —  D^G^u 1  =  IDGu l, 

«  ,  IDGunh)  =  «  -  xIDGwl ,  IDG{wl  -  xIDGw D) 

=  «,  IDGwl)  +  x2  (IDGwl ,  (. IDGfwl )  -  X  «  -  {IDGfwh) 

-  X  {IDGwnh  ,  IDGwnh) 

=  {< ,  IDGwl)  +  X2  {IDGwl ,  IDGwl)  ~  X2  {IDGwl ,  {IDG)DhGhwD 

-  2X  {IDGwl,  IDGwl) 

=  « ,  IDGwl)  -  X(2  -  X)  {IDGwl ,  IDGwl)  ~  X2  {IDGwl ,  {IDG)DhGhwnh)  , 
which  establishes  (5.45). 

To  establish  (5.47)  we  begin  with  (5.29).  Taking  the  inner-product  of  both  sides  with  respect  to  IDGe , 
(e,  IDGe)  =  (e,  IDGe)  -  x{IDGe,  IDGe)  +  x{IDGun,IDGe) , 

i.e. 

2  llell/.DG  +  2  ll^ll 7z?g  —  2  II e  £ll ?Z3G  =  Ikll/DG  ~  x{IDGe,IDGe)  +  x{IDGun ,IDGe) . 

Thus, 

IMI/dg  =  IMI/dg  ~  lle~  £\\2idg  +  2 x{IDGe,  IDGe)  —  2 x{IDGun ,  IDGe) .  (5.48) 
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In  addition,  rearranging  (5.29)  we  have 


e  —  e  =  —xIDGe  +  xIDGun 

and  thus, 

lle  -  42idg  =  ((e  -  £).  IDG(e  -  e))  (5.49) 

=  ((- XIDGe  +  XIDGun ) ,  IDG(-XIDG£  +  XIDGun )) 

=  X2\\IDGs\\2idg  +  x2\\IDGun\\2IDG  -  2X2(IDGs  ,  ( IDG)2un ) .  (5.50) 

Substituting  (5.50)  into  (5.48)  and  rearranging 

II£IIjdg  =  IMI/dg  ~  X2\\IDGe\\2IDG  —  x2 1 1 IDGun \ \ 2IDG  +  2 x2(IDGe ,  ( IDG)2un ) 

+  2 x(IDGe,IDGe)  -  2x{IDGun ,  IDGe) .  (5.51) 

Note  that 

2 x(IDGe,IDGe)  -  X2\\IDGe\\2DG  =  2X(IDGe,  IDGs)  -  X2(IDGe,(IDG)(IDG)e) 

=  (2X  -  x2){IDGe ,  IDGe)  +  X2(IDGe,  DhGh(IDG)e) ,  (5.52) 

2 x2(IDGe,  ( IDG)2un )  <  ^x2 (IDGe,  IDGe)  +  2 x2((IDG)2un  ,  {IDG)2un) ,  (5.53) 

2X(IDGun,  IDGe)  <  ^X(IDGe,IDGe)  +  2X(IDGun,IDGun) .  (5.54) 

Thus,  using  (5.52)-(5.54)  in  (5.51),  we  obtain 

\\42IDG  >  INI 2DG  +  |x(l  -  X)\\IDGe\\2  +  X2(IDGe,  DhGh(IDG)e) 

2X(1  +  X)\\IDGun\\2  -  x2\\IDGun\\2IDG. 

For  the  final  result  we  begin  again  with  e  —  e  =  — x/UGe  +  x/HGu".  Take  the  inner  product  with  e, 
use  the  polarization  identity  on  the  (e,  e)  term,  multiply  by  2  and  simplify.  This  gives 

||e||2  +  2 X{IDGe,e)  =  \\e\\2  +  \\e  -  e\\2 +  2X(IDGu,e).  (5.55) 

Taking  norms  we  also  have  ||e  —  e\\2  =  \\xIDGe  —  xIDGun\\2  =  y2(/I>G£  —  IDGun,IDGe  —  IDGun). 
Inserting  this  into  the  RHS  of  (5.55),  expanding  and  collecting  terms  yields 

llel|2+ 

x{(2  -  x)  ( IDGe,e )  +  x(IDGe,DGe)  -  x\\IDGu\\2  -  2(1  DGu,  e)  +  2 x(IDGe,  IDGu)} 

=  IN|2- 

For  the  last  two  terms  inside  the  braces  we  have,  using  operator  weighted  Cauchy-Schwarz  inequalities, 

2{IDGu,e)  <  2(1  DGu, u)  +  ^( IDGe,e ), 

2 X(IDGe,  IDGu)  =  2X(IDGe,  u)  -  2x(IDGe,  DGu) 

<  2 x(IDGu,u)  +  ~x{IDGe,e)  +  x(IDGe,  DGe)  +  x(IDGu,  DGu). 

Inserting  these  two  into  the  braces  and  collecting  terms,  we  have,  as  claimed: 

||e||2  <  ||e||2  -  ^x(l  -  x)(IDGe,e) 

+  [x2\\IDGu\\2  +  2x(l  +  x)\(I DGu, u)  +  x2 {IDGu,  DGu)\  . 

□ 
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6.  Numerical  Experiments.  In  this  section,  we  present  numerical  experiments  to  test  the  algorithms 
presented  herein.  Using  the  Green-Taylor  vortex  problem,  we  confirm  the  predicted  convergence  rates  of 
the  previous  section,  and  also  use  it  to  compare  the  accuracy  of  the  algorithms.  Further  testing  is  then 
performed  using  the  flow  around  a  cylinder  benchmark  problem.  We  use  the  software  FreeFEM++  [HePi] 
to  run  the  numerical  tests  which  provided  a  multi-frontal  Gauss  LU  factorization  for  the  linear  solver.  The 
scheme  presented  in  Algorithm  4.1  is  discretized  in  space  using  the  finite  element  method  with  Taylor- 
Hood  elements  (continuous  piecewise  quadratic  polynomials  for  the  velocity  and  continuous  linears  for  the 
pressure).  The  nonlinear  system  was  solved  by  a  fixed  point  iteration.  For  the  Stokes  filter,  used  in  all  the 
computations,  we  applied  the  same  boundary  conditions  as  given  for  the  problem  being  solved. 


m 

||  |^  ^ h  1 1|  oo,0 

rate 

|||Vlt  -  Vlt/,  2,0 

rate 

16 

2.55117- 10"2 

3.36582 

24 

2.05028  •  10"2 

0.54 

2.27963 

0.46 

32 

1.35885  •  10-2 

1.43 

1.6152 

1.18 

40 

8.1897-  KT3 

2.27 

1.18214 

1.73 

48 

5.01347  •  10-a 

2.69 

8.85739  •  10_i 

2.01 

56 

3.20674  •  10"a 

2.90 

6.76137-  10"1 

2.20 

64 

2.14058  •  10-a 

3.03 

5.24409  •  10"1 

2.33 

72 

1.48243  •  10"a 

3.12 

4.12513-  10"1 

2.44 

Table  6.1 

Errors  and  convergence  rates  for  order  of  deconvolution  N  =  I . 


m 

||  |^  ^ h  1 1|  oo,0 

rate 

1  Vw  — 

2,0 

rate 

16 

2.72906  •  10"2 

3.09567 

•  10_i 

24 

2.64931  •  10"2 

0.07 

3.01949 

•  10"1 

0.06 

32 

2.5414-  KT2 

0.14 

2.91435 

•  10"1 

0.12 

40 

2.40916- 10"2 

0.24 

2.78326 

•  io-1 

0.21 

48 

2.2577  •  KT2 

0.36 

2.63147 

•  10"1 

0.31 

56 

2.09387  •  10-2 

0.49 

2.46598 

•  io-1 

0.42 

64 

1.9251  •  KT2 

0.63 

2.29443 

•  10"1 

0.54 

Table  6.2 

Errors  and  convergence  rates  for  order  of  deconvolution  N  =  0. 


m 

||  l^1  ^ h  1 1|  oo,0 

rate 

III  Vlt  -  Vlift  ||  2,0 

rate 

16 

2.54983 

•  10"2 

2.93004  •  10"1 

24 

2.0471  • 

IO"2 

0.54 

2.42837  •  10”1 

0.46 

32 

1.35493 

•  10"2 

1.43 

1.72972  •  10"1 

1.18 

40 

2.40916 

•  10"3 

2.27 

1.17613-  10"1 

1.73 

48 

4.99331 

•  10"3 

2.69 

8.1466  •  10"2 

2.01 

56 

3.19319 

•  10"3 

2.90 

5.80612  •  10-2 

2.18 

64 

2.13123 

•  10"3 

3.03 

4.25236  •  10"2 

2.33 

Table  6.3 

Errors  and  convergence  rates  for  order  of  deconvolution  N  =  1  with  relaxation  \  =  At. 


6.1.  Convergence  Rate  Verification.  Our  first  test  is  designed  to  test  (and  does  confirm)  the  pre¬ 
dicted  rates  of  convergence.  The  problem  of  simulating  decay  of  the  Green-Taylor  vortex,  [GT37],  [T23],  is 
an  interesting  test  problem  in  which  the  true  solution  is  known.  It  was  used  as  a  numerical  test  in  Chorin 
[Cho68],  Tafti  [Tafti]  and  John  and  Layton  [JL02].  For  a  very  insightful  and  detailed  analysis  of  the  problem 
for  LES  models  see  Barbato,  Berselli  and  Grisanti  [BBG07]  and  Berselli  [B05].  The  prescribed  solution  in 
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SI  =  (0, 1)  x  (0, 1)  is  given  by 

Ui(x,y,t)  =  —  cos(w7ra;)  sin(w7ri/)e_2“  77  t^T  , 

r)  2  2./ 

U2(x,y,t)  =  sin(w7rx)  cos(w7rj/)e“  “  77  /r  , 
p(x,y,t )  =  —  ^-(cos(2w7ra:)  +  cos(2w7ry))e_2a’  77  t/T  . 

When  the  relaxation  time  r  =  i?e,  this  is  a  solution  of  the  NSE  with  /  =  0,  consisting  of  an  oj  x  u>  array  of 
oppositely  signed  vortices  that  decay  as  t  — »  oo. 

In  our  test,  with  Dirichlet  boundary  conditions  for  the  velocity,  we  choose  w  =  1,  At  =  0.005,  T  =  1, 
Reynolds  number  i?e  =  100  and  6  =  h  =  1/m,  where  m  is  the  number  of  subdivisions  of  the  interval  (0, 1). 
The  results  for  the  approximation  method  described  in  Algorithm  4.1  are  presented  in  Table  6.1,  using  order 
of  deconvolution  N  =  1  without  relaxation  (i.e.  X  =  0)  and  in  Table  6.3  for  N  =  1  with  relaxation  for 
X  =  At.  Results  using  the  simple  averaging  filter,  i.e.  deconvolution  with  order  N  =  0  ,  are  presented 
in  Table  6.2.  The  convergence  rate  is  calculated  from  the  error  at  two  successive  values  of  h  in  the  usual 
manner  by  postulating  e(h)  =  Ch 13  and  solving  for  /3  via  (3  =  ln(e(h\)  /  e(hf))  / ln(/ii/ti2). 

From  the  tables  we  see  the  convergence  rate  approaches  the  second  order  predicted  for  ||  |Vu  —  Vu/j|  ||2,o 
and  we  also  see  what  appears  to  be  an  L2  lift  for  ||  |it  —  Uh\  ||oo,o  for  order  of  deconvolution  N  =  1.  The  method 
with  order  of  deconvolution  N  =  0  ,  has  much  larger  errors  and  slower  rates  of  convergence,  as  expected. 
From  this  test  it  is  clear  that  deconvolution  makes  an  important  contribution  to  improving  the  accuracy  of 
the  approximation. 

6.2.  Flow  around  a  cylinder.  Our  next  numerical  illustration  is  for  two  dimensional  under-resolved 
flow  around  a  cylinder.  We  compute  values  for  the  maximal  drag  Cd,max  and  lift  ci^max  coefficient  at  the 
cylinder,  and  for  the  pressure  difference  A p[t)  between  the  front  and  back  of  the  cylinder  at  the  final  time 
T  =  8.  This  is  a  well  known  benchmark  problem  taken  from  Shafer  and  Turek  [ST96]  and  John  [ J 04] .  It 
is  not  turbulent  but  does  have  interesting  features.  The  flow  patterns  are  driven  by  the  interaction  of  a 
fluid  with  a  wall  which  is  an  important  scenario  for  industrial  flows.  This  simple  flow  can  be  difficult  to 
simulate  successfully  by  a  model  with  sufficient  regularization  to  handle  higher  Reynolds  number  problems. 
The  domain  is  presented  in  Figure  6.1. 

The  time  dependent  inflow  profile  is 

ui(0,  y,t)  =  Mi (2.2,  y,  t)  =  jj-|j^sin(7rt/8)t/(0.41  -  y) , 

U2(0,y,t)  =  112(2.2,  y,  t)  =  0. 

No  slip  boundary  conditions  are  prescribed  along  the  top  and  bottom  walls,  “do-nothing”  at  the  outflow  and 
the  initial  condition  is  u(x,y,0)  =  0.  The  viscosity  is  v  —  10-3  and  the  external  force  /  =  0.  The  Reynolds 
number  of  the  flow,  based  on  the  diameter  of  the  cylinder  and  on  the  mean  velocity  inflow  is  0  <  Re  <  100. 
A  mesh  with  62757  number  of  degrees  of  freedom  is  used  for  all  simulation  for  a  clear  comparison  of  the 
different  algorithms  presented  in  this  report.  The  filter  radius  is  chosen  as  the  length  of  the  cylinder  divided 
by  the  number  of  mesh  points  around  the  cylinder. 
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Fig.  6.2.  The  velocity  at  t  =  2,  4 ,  5,  6,  7,  and  8  of  Algorithm  4-1,  with  N  =  1  and  \  =  At  =  0.01. 


Fig.  6.3.  The  development  of  Cd(t),  ci(t)  and  A pit)  of  Algorithm  4-1  with  N  =  1  and  x  =  At  =  0.01. 


For  this  setting,  it  is  expected  that  as  the  flow  increases  from  time  t  =  2  to  t  =  4  two  vortices  start  to 
develop  behind  the  cylinder.  Between  t  =  4  and  t  =  5,  the  vortices  separate  from  the  cylinder,  so  that  a 
vortex  street  develops,  and  they  continue  to  be  visible  through  the  final  time  t  =  8.  This  can  be  seen  in 
Figure  6.2. 

The  evolutions  of  Cd,max,  citmax  and  A p  in  time  are  presented  in  Figure  6.3. 

For  the  computation  of  drag  and  lift  coefficients  we  used  the  one  dimensional  method  described  by  John 
in  [ J 04] .  Results  on  the  computations  of  maximal  drag  and  lift  coefficients  and  pressure  drop,  for  N  =  1, 


21 


are  presented  in  Table  6.4.  The  following  reference  intervals  are  given  in  [ST96]: 

Cdtax  e  [2-93,  2.97],  c\%ax  G  [0.47,  0.49],  A tf*  G  [-0.115,  -0.105] 


relax,  coeff. 

At 

t(Q2,matc) 

Cd,max 

t(Q,matc) 

Cl, max 

Ap(8s) 

0.04 

3.96 

2.87862 

6.12 

0.36572 

-0.103274 

X  =  o 

0.02 

3.94 

2.85123 

5.98 

0.420198 

-0.103153 

0.01 

3.94 

2.80214 

5.96 

0.409368 

-0.106276 

0.04 

3.96 

2.94233 

6.12 

0.392735 

-0.102161 

X  =  At 

0.02 

3.94 

2.94317 

5.98 

0.458007 

-0.106762 

0.01 

3.94 

2.94352 

5.93 

0.479286 

-0.110899 

Table  6.4 

Results  for  drag/lift  coefficients  and  pressure  difference  for  deconvolution  order  N  =  1. 

The  computations  in  Table  6.4  show  that  evolve- filter- deconvolve-relax,  Algorithm  4.1,  computes  the 
drag  and  lift  coefficients,  and  the  pressure  difference,  within  the  benchmark  intervals,  and  illustrates  the 
positive  role  of  using  relaxation  in  the  approximation  algorithm.  Computations  were  also  performed  for 
higher  Reynolds  number  corresponding  to  v  =  10~4.  For  this  case  the  direct  approximation  approach 
(no  regularization)  failed  (i.e.  the  fixed  point  iteration  for  the  discretization  of  the  nonlinear  term  of  the 
Navier-Stokes  equations  stopped  to  converge  around  T  =  6),  whereas  computations  using  Algorithm  4.1  were 
successful  for  N  =  1,  see  Figure  6.4.  There  are  no  benchmark  intervals  for  lift  and  drag  coefficients  for  this 
case.  Some  quality  of  the  approximation  can  be  assessed  by  the  appropriate  appearance  and  evolution  of 
the  vortex  street  in  the  simulation.  The  direct  approximation  without  regularization  is  much  more  sensitive 
to  the  instabilities  that  appear  in  the  flow  at  higher  Re  than  Algorithm  4.1,  which  produced  a  clear  vortex 
street  given  in  Figure  6.4  for  N  —  1. 

7.  Conclusions.  In  2002  R.  Peyret  wrote  in  the  fundamental  book  on  spectral  methods  for  the  NSE 
[P02]  that 

“...  filters  must  be  used  with  care  and  their  effect  evaluated  with  precision  when  tuning 
the  parameters... The  application  of  a  filter  in  the  course  of  time  integration  (especially  if  it 
is  often  applied)  may  be  dangerous.” 

We  have  found  that  simple  filtering  at  each  step  is  indeed  dangerous  as  it  introduces  large  amounts 
of  numerical  diffusion.  On  the  other  hand,  our  precise  analysis  has  also  shown  that  with  care,  filtering 
plus  deconvolution  plus  relaxation  stabilizes  marginally  resolved  scales  and  does  not  over  diffuse.  This 
combination  can  be  an  invaluable  tool,  with  many  algorithmic  advantages,  for  enhancing  stability  without 
degrading  accuracy  through  an  extra  uncoupled  and  modular  step. 

We  have  only  considered  the  simplest  filter  and  deconvolution  operators  that  fit  our  mathematical  tools. 
There  are  very  many  other  possible  filters  and  deconvolution  operators  that  can  be  tested  and  await  analysis. 
From  one  point  of  view,  the  process  Filter  — >  Deconvolve  simply  generates  another  filter  which  is  closer  to 
spectral  cutoff  and  which  can  be  applied  by  a  sequence  of  simpler  filter  steps.  Thus,  both  investigation  of 
other  filters  and  study  of  choices  of  relaxation  parameters  are  important  next  steps. 
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